Analysis parameters:

iGraph Parameters

connectom betas: refit full dataset using nested best lambda
connection counts (duplicates): YES

Load Data

Load intermediate data after running worksapce_ml.RData

load(file = "../data/__cache__/worksapce_ml.RData")

USE_NESTED <- TRUE 
SAVE_PLOT <- FALSE

Brain Network Analysis (Power)

A Matrix

To calculate the averaged corr matrix A

find the Fisher’s Z values of the corresponding Pearson correlation coefficients
Average them
Find the reverse Fisher’s Z transform of that average value.

C.z <- FisherZ(C)
C.zmean <- matrix(colMeans(C.z), nrow=264, ncol = 264)
A <- FisherZInv(C.zmean)
A.vec <- as.vector(A)

W Matrix

Calculate W from betas and A Matrix

connectom2matrix <- function(connectome, w) {
  empty_mat <- matrix(0, 264, 264, dimnames = list(paste0("X", 1:264), paste0("X", 1:264))) 
  empty_mat[connectome$index] <- connectome$W
  return(empty_mat)
}

# convert a 264*264 matrix back to connectom df
matrix2connectom <- function(mat, connectome, col_name) {
  connectome$temp = mat[connectome$index]
  connectome <- rename(connectome, !!col_name := temp)
  return(connectome)
}

# make the matrix symmetric
make_symmetric <- function(m) {
  # lower.tri is 0.0
  m[lower.tri(m)] <- t(m)[lower.tri(m)]
  return(m)
}

power_atals <- power2011 %>% 
  rename(ROI.Name = ROI, x.mni=X, y.mni=Y, z.mni=Z, network=NetworkName) %>% 
  mutate(ROI.Name=as.integer(ROI.Name), index = as.integer(ROI.Name),
         x.mni=as.integer(x.mni), y.mni=as.integer(y.mni), z.mni=as.numeric(z.mni)) %>%
  dplyr::select(ROI.Name, x.mni, y.mni, z.mni, network, index)

check_atlas(power_atals)

if(USE_NESTED) {
  connectome_data <- connectome_nested
} else {
  connectome_data <- connectome
}

Wconnectome <- connectome_data %>%
  mutate(A = A.vec[connectome_data$index], W = A*Beta)
  #separate(connection, c("connection1", "connection2"))%>%
  #separate(network, sep = "-", c("network1", "network2"), remove = F) %>%
  #filter(str_detect(network, pattern = "-1-")) %>%
         #network1 = ifelse(str_detect(network, pattern = "-1-"), -1, network1)) %>%
  #mutate(connection_type = ifelse(network1==network2, "Within", "Between"))

W_mat <- matrix(0, ncol = 264, nrow = 264)
W_mat[Wconnectome$index] <- Wconnectome$W
W_mat <- make_symmetric(W_mat) #CHECKED correct W_mat


rownames(W_mat) = power_atals$network
colnames(W_mat) = power_atals$network

Declarative Network

power_color = power2011 %>% 
  filter(NetworkName!="Sensory/somatomotor Mouth") %>%
  select(NetworkName, Color) %>% 
  distinct() 

colors <- c(Uncertain = power_color$Color[9], #power_color$Color[1], 
            `Sensory/somatomotor Hand` = power_color$Color[6], #power_color$Color[2], 
            `Cingulo-opercular Task Control` = power_color$Color[3],
            Auditory = power_color$Color[4],
            `Default mode` = power_color$Color[5], 
            `Memory retrieval?` = power_color$Color[2], #power_color$Color[6],
            `Ventral attention` = power_color$Color[7], 
            Visual = power_color$Color[8],
            `Fronto-parietal Task Control` = power_color$Color[1], #power_color$Color[9],
            Salience = power_color$Color[10],
            Subcortical = power_color$Color[11], 
            Cerebellar = power_color$Color[12], 
            `Dorsal attention` = power_color$Color[13])

colors <- c(Uncertain = power_color$Color[1], 
            `Sensory/somatomotor Hand` = power_color$Color[2],  
            `Cingulo-opercular Task Control` = power_color$Color[3],
            Auditory = power_color$Color[4],
            `Default mode` = power_color$Color[5], 
            `Memory retrieval?` = power_color$Color[6], 
            `Ventral attention` = power_color$Color[7], 
            Visual = power_color$Color[8],
            `Fronto-parietal Task Control` = power_color$Color[9], 
            Salience = power_color$Color[10],
            Subcortical = power_color$Color[11], 
            Cerebellar = power_color$Color[12], 
            `Dorsal attention` = power_color$Color[13])

if (SAVE_PLOT) {
  png(filename = "analysis_lasso_re_3connectom/figure-html/network_figure1.png",  bg = "transparent", 
    width = 800 , height = 800, 
    units = "px", res = 150)
  
  circos.clear()
  chordDiagram(W_mat, directional = FALSE, transparency = 0.5, self.link = 1, 
             symmetric = TRUE, scale = TRUE, reduce = FALSE, 
             annotationTrackHeight = mm_h(c(15, 1)),
             annotationTrack = c("grid"),
             grid.col = colors, col = ifelse(W_mat>0, "tomato", "#00000000")) 
  title("Predictive Group: Declarative (W > 0)", outer = F, cex.main = 1.5)
} else {
  circos.clear()
  chordDiagram(W_mat, directional = FALSE, transparency = 0.5, self.link = 1, 
             symmetric = TRUE, scale = TRUE, reduce = FALSE, 
             #annotationTrackHeight = mm_h(c(15, 1)),
             #annotationTrack = c("grid"),
             grid.col = colors, col = ifelse(W_mat>0, "tomato", "#00000000")) 
  title("Predictive Group: Declarative (W > 0)", outer = F, cex.main = 1.5)
  
  }

Procedural Network

if (SAVE_PLOT) {
  png(filename = "analysis_lasso_re_3connectom/figure-html/network_figure2.png",  bg = "transparent", 
    width =800 , height = 800, 
    units = "px", res = 150)

  circos.clear()
  chordDiagram(W_mat, directional = FALSE, transparency = 0.5, self.link = 1, 
               symmetric = TRUE, scale = TRUE, reduce = FALSE, 
               annotationTrackHeight = mm_h(c(15, 1)),
               annotationTrack = c("grid"),
               grid.col = colors, col = ifelse(W_mat<0, "steelblue", "#00000000"))
  title("Predictive Group: Procedural (W < 0)", outer = F, cex.main = 1.5) 
  } else {
  circos.clear()
  chordDiagram(W_mat, directional = FALSE, transparency = 0.5, self.link = 1, 
               symmetric = TRUE, scale = TRUE, reduce = FALSE, 
               #annotationTrackHeight = mm_h(c(15, 1)),
               #annotationTrack = c("grid"),
               grid.col = colors, col = ifelse(W_mat<0, "steelblue", "#00000000"))
  title("Predictive Group: Procedural (W < 0)", outer = F, cex.main = 1.5)
}

circos.clear()
chordDiagram(W_mat, directional = FALSE, transparency = 0.5, self.link = 1, 
             symmetric = TRUE, scale = TRUE, reduce = FALSE, 
             order = unique(power2011$NetworkName),  
             preAllocateTracks = 1,
             annotationTrackHeight = mm_h(c(10, 10)),
             annotationTrack = c("grid"),
             grid.col = colors, col = ifelse(W_mat>0, "tomato", "steelblue"))
title("Declarative + Procedural Network",  cex.main = 2)
circos.trackPlotRegion(track.index = 1, panel.fun = function(x, y) {
  xlim = get.cell.meta.data("xlim")
  ylim = get.cell.meta.data("ylim")
  sector.name = get.cell.meta.data("sector.index")
  circos.text(mean(xlim), ylim[1] + .1, sector.name, facing = "clockwise", niceFacing = TRUE, adj = c(0, 0.5))
  circos.axis(h = "top", labels.cex = 0.5, major.tick.percentage = 0.2, sector.index = sector.name, track.index = 2)
}, bg.border = NA)

Network and node desciptives

Next, we will look at Graph properties of two networks

###Graph Density

The proportion of present edges from all possible edges in the network.

# select cols
roi_links <- Wconnectome %>% dplyr::select(connection1, connection2, W, connection_type, network, network_names)
# rename cols
colnames(roi_links) <- c("from", "to", "weight", "connection_type", "network", "network_names")

roi_nodes <- power2011 %>% rename(id = ROI) %>%
  mutate(NetworkName=factor(NetworkName), 
         Color=factor(Color))
levels(roi_nodes$Color) <- sample(colors(T), 14) 

# create a graph
net <- graph_from_data_frame(d=roi_links, vertices=roi_nodes, directed=F) 
#g <- graph_from_adjacency_matrix(W_mat,, mode = "upper")

net.d <- net - E(net)[E(net)$weight<0]
net.p <- net - E(net)[E(net)$weight>0]

df.density <- data.frame("edge_density"=c(edge_density(net.d, loops=F), 
                            edge_density(net.p, loops=F),
                            edge_density(net, loops=F)), 
           "network" = factor(c("Declarative", "Procedural", "Full"), levels = c("Declarative", "Procedural", "Full")))

df.density %>% ggbarplot(x="network", y="edge_density", color = "white", fill=c("tomato", "steelblue", "gray"), width = 0.5, label = T, lab.nb.digits = 5) +
  theme_pander() + 
  ggtitle("Degree Edge Density")

Graph: Diameter

A network diameter is the longest geodesic distance (length of the shortest path between two nodes) in the network. In igraph, diameter() returns the distance, while get_diameter() returns the nodes along the first found path of that distance.

# make negative weights to positive
net.p.abs <- net.p
E(net.p.abs)$weight <- E(net.p.abs)$weight * (-1)

# make negative weights to positive
net.abs <- net
E(net.abs)$weight[E(net.abs)$weight<0] <- E(net.abs)$weight[E(net.abs)$weight<0] * (-1)


df.diameter <- data.frame("diameter"=c(diameter(net.d, directed=F), 
                        diameter(net.p.abs, directed=F),
                        diameter(net.abs, directed=T)), 
           "network" = factor(c("Declarative", "Procedural", "Full"), levels = c("Declarative", "Procedural", "Full")))

df.diameter %>% ggbarplot(x="network", y="diameter", color = "white", fill=c("tomato", "steelblue", "gray"), width = 0.5, label = T, lab.nb.digits = 5) +
  theme_pander() + 
  ggtitle("Network Diameter")

Graph: Centrality Degree

Centrality functions (vertex level) and centralization functions (graph level). The centralization functions return res - vertex centrality, centralization, and theoretical_max - maximum centralization score for a graph of that size. The centrality function can run on a subset of nodes (set with the vids parameter). This is helpful for large graphs where calculating all centralities may be a resource-intensive and time-consuming task.

Centrality is a general term that relates to measures of a node’s position in the network. There are many such centrality measures, and it can be a daunting task to wade through all of the different ways to measure a node’s importance in the network. Here, we will introduce just a few examples.

df.centrality <- data.frame("centr_degree"=c(centr_degree(net.d, normalized=T)$centralization, 
                        centr_degree(net.p, normalized=T)$centralization,
                        centr_degree(net, normalized=T)$centralization), 
           "network" = factor(c("Declarative", "Procedural", "Full"), levels = c("Declarative", "Procedural", "Full")))

df.centrality %>% ggbarplot(x="network", y="centr_degree", color = "white", fill=c("tomato", "steelblue", "gray"), width = 0.5, label = T, lab.nb.digits = 5) +
  theme_pander() + 
  ggtitle("Network Degree Centrality ")

Graph: Betweeness (Closeness)

Let’s now do the same for betweenness centrality, which is defined as the number of geodesic paths (shortest paths) that go through a given node. Nodes with high betweenness might be influential in a network if, for example, they capture the most amount of information flowing through the network because the information tends to flow through them. Here, we use the normalized version of betweenness.

Closeness (centrality based on distance to others in the graph) Inverse of the node’s average geodesic distance to others in the network.

df.sloseness <- data.frame("centr_clo"=c(centr_clo(net.d)$centralization, 
                        centr_clo(net.p)$centralization,
                        centr_clo(net)$centralization), 
           "network" = factor(c("Declarative", "Procedural", "Full"), levels = c("Declarative", "Procedural", "Full")))

df.sloseness %>% ggbarplot(x="network", y="centr_clo", color = "white", fill=c("tomato", "steelblue", "gray"), width = 0.5, label = T, lab.nb.digits = 5) +
  theme_pander() + 
  ggtitle("Network Closeness")

Graph: Distances

df.distance <- data.frame("mean_distance"=c(mean_distance(net.d, directed=F), 
                        mean_distance(net.p, directed=F),
                        mean_distance(net, directed=F)), 
           "network" = factor(c("Declarative", "Procedural", "Full"), levels = c("Declarative", "Procedural", "Full")))

df.distance %>% ggbarplot(x="network", y="mean_distance", color = "white", fill=c("tomato", "steelblue", "gray"), width = 0.5, label = T, lab.nb.digits = 5) +
  theme_pander() + 
  ggtitle("Network Closeness")

Graph: Assortativity

df.assortativity<- data.frame("assortativity_degree"=c(assortativity_degree(net.d, directed=F), 
                        assortativity_degree(net.p, directed=F),
                        assortativity_degree(net, directed=F)), 
           "network" = factor(c("Declarative", "Procedural", "Full"), levels = c("Declarative", "Procedural", "Full")))

df.assortativity %>% ggbarplot(x="network", y="assortativity_degree", color = "white", fill=c("tomato", "steelblue", "gray"), width = 0.5, label = T, lab.nb.digits = 5) +
  theme_pander() + 
  ggtitle("Network Assortativity Degree")

> Combined Table

options(scipen=999)
df.density %>% 
  select(network, edge_density) %>%
  left_join(df.diameter) %>%
  left_join(df.centrality) %>%
  left_join(df.sloseness) %>%
  left_join(df.distance) %>%
  left_join(df.assortativity) %>%
  rename("centrality_degree"=centr_degree, 
         "centrality_closeness"=centr_clo) %>%
  kable(format.args = list(scientific = TRUE, big.mark = ",", digit=3))

network	edge_density	diameter	centrality_degree	centrality_closeness	mean_distance	assortativity_degree
Declarative	9.22e-04	1.12e+00	6.68e-03	5.00e-05	1.14e+00	-1.85e-01
Procedural	1.56e-03	2.73e+00	1.37e-02	1.83e-04	1.83e+00	-8.19e-02
Full	2.48e-03	2.73e+00	1.27e-02	3.04e-04	2.56e+00	4.60e-02

Graph: Similarity Measurement

W values in Brain Connectom

#colors <- factor(power_atals$network)
#levels(colors) <- colors14
#power_atals$colors <- as.character(temp)
check_atlas(power_atals)
x1 <- W_mat
x1[x1<0] <- 0


p1 <- brainconn(atlas=power_atals, conmat=x1, node.color = power2011$Color, view = "top",
          node.size = igraph::degree(net.d)*2.5, all.nodes = TRUE, 
          edge.color = "tomato",  edge.color.weighted = FALSE, scale.edge.width = c(1,3), edge.alpha = 0.6,
          label.edge.weight = F,  show.legend = F,
          background.alpha = .3) #+ ggtitle("Strategy Predictability: W")
p2 <- brainconn(atlas=power_atals, conmat=x1, node.color = power2011$Color, view = "left",
          node.size = igraph::degree(net.d)*2.5, all.nodes = TRUE, 
          edge.color = "tomato",  edge.color.weighted = FALSE, scale.edge.width = c(1,3), edge.alpha = 0.6,
          label.edge.weight = F,  show.legend = F,
          background.alpha = .3) #+ ggtitle("Strategy Predictability: W")
p3 <- brainconn(atlas=power_atals, conmat=x1, node.color = power2011$Color, view = "back",
          node.size = igraph::degree(net.d)*2.5, all.nodes = TRUE, 
          edge.color = "tomato",  edge.color.weighted = FALSE, scale.edge.width = c(1,3), edge.alpha = 0.6,
          label.edge.weight = F,  show.legend = F,
          background.alpha = .3) #+ ggtitle("Strategy Predictability: W")

x2 <- W_mat
x2[x2>0] <- 0

p4 <- brainconn(atlas=power_atals, conmat=x2*-10, node.color = power2011$Color, view = "top",
          node.size = igraph::degree(net.p)*2.5, all.nodes = TRUE, 
          edge.color = "steelblue",  edge.color.weighted = FALSE, 
          scale.edge.width = c(1,3), edge.alpha = 0.6,
          label.edge.weight = F,  show.legend = F,
          background.alpha = .3)  



p5 <- brainconn(atlas=power_atals, conmat=x2*-10, node.color = power2011$Color, view = "left",
          node.size = igraph::degree(net.p)*2.5, all.nodes = TRUE, 
          edge.color = "steelblue",  edge.color.weighted = FALSE, 
          scale.edge.width = c(1,3), edge.alpha = 0.6,
          label.edge.weight = F,  show.legend = F,
          background.alpha = .3) 

p6 <- brainconn(atlas=power_atals, conmat=x2*-10, node.color = power2011$Color, view = "back",
          node.size = igraph::degree(net.p)*2.5, all.nodes = TRUE, 
          edge.color = "steelblue",  edge.color.weighted = FALSE, 
          scale.edge.width = c(1,3), edge.alpha = 0.6,
          label.edge.weight = F,  show.legend = F,
          background.alpha = .3)

if (SAVE_PLOT) {
  png(filename = "analysis_lasso_re_3connectom/figure-html/network_figure4.png",  bg = "transparent", 
      width =500, height = 500, units = "px", res = 150) }
p1

if (SAVE_PLOT) {
  png(filename = "analysis_lasso_re_3connectom/figure-html/network_figure5.png",  bg = "transparent", 
      width =500, height = 500, units = "px", res = 150)} 
p2

if (SAVE_PLOT) {
  png(filename = "analysis_lasso_re_3connectom/figure-html/network_figure6.png",  bg = "transparent", 
      width =500, height = 500, units = "px", res = 150)}
p3

if (SAVE_PLOT) {
  png(filename = "analysis_lasso_re_3connectom/figure-html/network_figure7.png", 
      bg = "transparent", width =500, height = 500, units = "px", res = 150)}
p4

if (SAVE_PLOT) {
  png(filename = "analysis_lasso_re_3connectom/figure-html/network_figure8.png",  bg = "transparent", 
      width =500, height = 500, units = "px", res = 150)}

p5

if (SAVE_PLOT) {png(filename = "analysis_lasso_re_3connectom/figure-html/network_figure9.png",  bg = "transparent", 
                    width =500, height = 500, units = "px", res = 150)}

p6

Distribution of connections

Look at the distribution of network in two groups

DUPLICATE <- TRUE

c1 <- Wconnectome %>% filter(W>0) %>% mutate(roi = as.integer(connection1)) %>% dplyr::select(roi) %>%  unlist()
c2 <- Wconnectome %>% filter(W>0) %>% mutate(roi = as.integer(connection2)) %>% dplyr::select(roi) %>%  unlist()


c3 <- Wconnectome %>% filter(W<0) %>% mutate(roi = as.integer(connection1)) %>% dplyr::select(roi) %>%  unlist()
c4 <- Wconnectome %>% filter(W<0) %>% mutate(roi = as.integer(connection2)) %>% dplyr::select(roi) %>%  unlist()


if (DUPLICATE) {
  c12 <- c(c1, c2)
  c34 <- c(c3, c4)
} else {
  c12 <- unique(c(c1, c2))
  c34 <- unique(c(c3, c4))
}

df.c1 <- power2011[c12,] %>%
  mutate(NetworkName = factor(NetworkName)) %>%
  count(NetworkName, name = "count", .drop = F)%>%
  right_join(power2011 %>% dplyr::select(NetworkName, Color) %>% distinct(), on="NetworkName") %>%
  mutate(count=as.integer(ifelse(is.na(count), 0, count))) %>%
  arrange(NetworkName)

p7 <- power_color %>% 
  left_join(df.c1) %>%
  ggplot(aes(x=count, y=NetworkName)) +
  geom_col(fill = power_color$Color) +
  scale_fill_manual(guide = guide_legend(reverse = F), 
                    values = df.c1$Color, 
                    labels =  df.c1$Color) +
  scale_x_reverse() +
  theme_pander() +
  ggtitle("Distribution of connections", subtitle = "Declarative Network") +
  theme(legend.position = "right", 
        axis.text.y = element_blank(), 
        axis.title.y = element_blank(),
        plot.title = element_text(size = 20),
        axis.text = element_text(size = 20),
        legend.title = element_text(size = 20),
        legend.text = element_text(size = 20))

p7

# count number of networks included in
df.c2 <- power2011[c34,] %>%
  mutate(NetworkName = factor(NetworkName)) %>%
  count(NetworkName, name = "count", .drop = F) %>%
  right_join(power2011 %>% dplyr::select(NetworkName, Color) %>% distinct(), on="NetworkName") %>%
  mutate(count=ifelse(is.na(count), 0, count)) %>%
  arrange(NetworkName) 



p8 <- power_color %>% 
  left_join(df.c2) %>%
  ggplot(aes(y = count, x=NetworkName)) +
  geom_col(fill = power_color$Color) +
  coord_flip() +
  scale_fill_manual(guide = guide_legend(reverse = F), values = power_color$Color) +
  theme_pander() +
  ggtitle("Distribution of connections", subtitle = "Procedural Network") +
  theme(legend.position = "right", 
        axis.text.y = element_blank(), 
        axis.title.y = element_blank(),
        plot.title = element_text(size = 20),
        axis.text = element_text(size = 20),
        legend.title = element_text(size = 20),
        legend.text = element_text(size = 20))



p8

if (SAVE_PLOT) {
  png(filename ="analysis_lasso_re_3connectom/figure-html/network_figure10.png",  bg = "transparent", 
    width = 2000, height = 2000, 
    units = "px", res = 300)}

ggbarplot(df.c2 , x="NetworkName", y="count", fill = "NetworkName", color="white",
          palette = df.c2$Color, #order = "count",    
          legend = "right",
          rotate = TRUE, ggtheme = theme_pander(),
          title = "Distribution of connections") +
labs(fill = "Network Name")

# ggbarplot(df.c2 , x="NetworkName", y="count", fill = "NetworkName", color="white",
#           palette = df.c2$Color, #order = "count",    
#           
#           rotate = TRUE, ggtheme = theme_pander(),
#           title = "Distribution of connections") +
#   scale_y_continuous(limits = c(0,15), breaks=c(0,5, 10, 15)) +
#   scale_fill_discrete(breaks=rev(df.c2$NetworkName), ) +
#   scale_fill_discrete(guide = guide_legend(reverse=TRUE)) +
#   theme(legend.position = "right",
#         axis.text.y = element_blank(),
#         axis.title.y = element_blank(),
#         plot.title = element_text(size = 20),
#         axis.text = element_text(size = 20),
#         legend.title = element_text(size = 20),
#         legend.text = element_text(size = 20))

if (SAVE_PLOT) {
  png(filename = "analysis_lasso_re_3connectom/figure-html/network_figure10.png",  bg = "transparent", 
    width =3000, height = 1000, 
    units = "px", res = 150)}

ggarrange(p7, NULL, NULL, NULL, p8, 
          #labels = c("A", "B", "C"),
          ncol = 5, nrow = 1, align = "h", 
          common.legend = TRUE, legend = "bottom", widths = c(1,.8,.8,1))

Chi-sq Test

chisq.test(df.c1$count, df.c2$count, simulate.p.value = T, p = noi_stats$N)

## 
##  Pearson's Chi-squared test with simulated p-value (based on 2000
##  replicates)
## 
## data:  df.c1$count and df.c2$count
## X-squared = 98, df = NA, p-value = 0.2994

Analysis Lasso-Re (Connectom)

Cher Yang

2023-02-13